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Abstract. - We argue that a system of straight rigid rods of length k on square lattice with 
only hard-core interactions shows two phase transitions as a function of density p for k > 7. 
The system undergoes a phase transition from the low-density disordered phase to a nematic 
phase as p is increased from at p = pd, and then again undergoes a reentrant phase transition 
from the nematic phase to a disordered phase at p = pc2 < 1- 



Introduction. - The study of systems of long rod-like molecules in solution with only 
excluded volume interaction has a long history. Onsager had noted that such a system would 
show long-range orientational order at large enough density [1] . Flory studied a lattice model 
of long rod-like molecules, and based on a mean-field approximation, argued that the lattice 
model would also show an isotropic- nematic phase transition as a function of density [2]. 
Zwanzig studied a model of hard rods in the continuum, but with only a finite number of 
allowed orientations [3]. For the continuum problem, there is general agreement that in three 
dimensions, the isotropic-nematic transition occurs at large enough density, if the length-to- 
width ratio of the molecules is large enough [4]. In the case of two-dimensional systems of 
hard needles in a plane, there is no spontaneous breaking of continuous symmetry, but the 
high-density phase has orientational correlations that decay as a power- law [5,6]. 

The situation is much less clear in the case of lattice models of straight hard rods of length 
k. The only analytically soluble case is that of dimers (k — 2), for which it is known that the 
orientational correlations have an exponential decay for all densities, except the full-coverage 
limit, which shows power-law correlations in all dimensions [7,8]. Power-law correlations at 
p = 1 were also seen in numerical study on the square lattice for fc = 3 [9], and k = 4 [10]. 
Early Monte Carlo studies of semi-flexible fc-mers on a square lattice showed no long-range 
orientational order for any density or temperature [11]. DeGennes and Prost, in 1995, in their 
book on liquid crystals, wrote "it is not quite certain that such a lattice model, even for large 
p (long molecule), will lead to a transition" [12]. While there has been some work on this 
general problem since then, the question still remains unsettled [13] . 

In this paper, we present the results of our Monte Carlo study of this model on the 
square lattice. We find very strong numerical evidence that the system shows nematic order 
at intermediate densities for k > 7. We are not able to study systems near p = 1 in our 
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Fig. 1 - Schematic representation of different phases of the hard rod system on a square lattice, 
(a) The low-density disordered phase (b) Intermediate density ordered phase (c) the high density 
disordered phase. 



simulations, because the relaxation time increases very fast with density. However, there is 
fairly convincing numerical evidence (but no proof), that there is no long-range order at p = 1, 
and so there must be a second phase transition as density is increased. To the best of our 
knowledge, this second phase-transition has not been discussed in published literature before. 
We estimate the critical density for the second transition. The different phases are shown 
schematically in Fig. 1. 

Model. - Our model consists of a system of hard straight rods, each of which occupies k 
consecutive sites on a lattice. A rod has only two possible orientations, horizontal or vertical. 
We will call each such rod a fc-mer. The only interaction between different rods is hard core 
exclusion: no site can be occupied by more than one fc-mer. 

We take the lattice to be a two-dimensional square lattice of size L x L. Let N{nh,nv) be 
the number of distinct configurations with nh horizontal and Uy vertical fc-mers. We define Zh 
and Zv as the activities for the horizontal and vertical fc-mers, and define the grand partition 
function 

n{zh,z,,)^ N{nH,ny)zl^z'^^. (1) 

The average values of rih and n„ in the grand canonical ensemble are obtained by taking 
partial derivatives of Vl with respect to Zh and z^. We define the mean density p as the 
fractional number of sites occupied by fc-mers. The nematic order parameter Q is defined as 



Q = Lirriz \z, LirriL^ao-r— Ti (2) 

where the limit is taken with Zy tending to Zh from above. We will usually be interested in 
the and denote the common value by z without a subscript. 

Results of the Numerical Simulations. - We have studied this problem by Monte Carlo 
simulations using a deposition-evaporation algorithm: At each time-step, we make an attempt 
to deposit a rod with probability p, or evaporate one of the existing rods with probability 1 —p. 
If it is a deposition attempt, we choose a horizontal or vertical orientation with probability 
1/2 each, and pick one of the sites of the lattice at random, and attempt to deposit the rod 
there in the chosen direction. The attempt is successful if that site, and the next fc — 1 sites in 
the chosen direction are all empty. Otherwise, the attempted move is rejected. We use period 
boundary conditions in both the horizontal and vertical directions. It is easy to see that this 
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Fig. 2 - Snapshots of time evolution of fc-mers, = 8, to the nematically ordered state for p = 
.8, starting from an empty lattice . Initially, the algorithm produces a disordered state, which 
slowly relaxes to the almost perfectly orientationally ordered state. The frames a-f are taken after 
1000, 4000, 8000, 12000, 20000 and 36000 MCS. The asymptotic value of the order parameter is g = 
0.989 and density, p = 0.758. 



deposition-evaporation dynamics is ergodic and satisfies detailed balance, with the parameter 
p related to activity z by 



PP 



(3) 



The deposition-evaporation dynamics does not conserve the number of rods, n and the 
number of vertical and horizontal rods n„ and Uh vary with time. In our Monte Carlo studies, 
we studied different values of k, from 2 to 12. We varied the deposition rate p and monitored 
the density p and the order parameter Q. 

Our main result is the evidence for existence of a orientationally ordered phase for in- 
termediate densities of rods for k > 7. The existence of the nematic phase is very clearly 
established in our simulations. A typical run is shown in Fig. [21 where we have taken k — 8, 
on a 80 X 80 lattice, for p — 0.8, starting from all empty initial state. In the beginning, both 
horizontal and vertical deposition events are equally likely, and the density builds up to nearly 
equilibrium value ~ .75, but the order parameter Q remains nearly zero. By the time 1000 
MCS, the density has nearly stabilized to its equilibrium value, but the order parameter is 
nearly zero. With time, the size of locally ordered domains grow, until one of the domains 
nearly fills the full lattice. For our lattice, this happens by around 16000 MCS. At this stage, 
the order parameter is nearly 1, and it shows only very small fiuctuations away from this value 
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Fig. 3 - Order parameter Q as a function of densities p is shown for different k and L. 

for the full length of simulation ( about lO'^ MCS). The average value of Q in this run was 
0.989. Of course, the ordering is equally likely to be horizontal or vertical. 

We have tested the results for different initial configurations, namely, all empty, all vertical, 
all horizontal, half filled, etc. to ensure sufficient time has been given for equilibration, and 
the averages calculated are independent of initial conditions. 

In Fig. Owe show that the order parameter Q as a function of density p for different length 
of rods fc = 6, 7,8, 10 and two system sizes L = 90 and 120. We did not find any isotropic- 
nematic transition for short rods fc < 6. However, as our simulations only cover the range 
< p < .85, we can only conclude that if kmin the lowest value of k that shows a nematic 
phase, then kmin < 7. 

For long rods fc > 7, a nematic phase is observed for intermediate densities. The critical 
density for the onset of nematic order decreases with increasing k. This is consistent with 
the theoretical expectation that for very large k, the critical density should be related to that 
found in the problem of ultra-thin needles in a continuum. If the critical density of needles 
of length 1 in continuum is Ac, the critical density for the isotropic to nematic transition pd 
should be approximately equal to Ac/k, for large k [Ac « 6, from simulations]. 

The point of onset of the nematic order is not so easy to determine from the data. In 
Fig. 3] we show the distribution of —n^ for fc = 10 and L — 120 for different densities p. For 
p < 0.36 the distribution is centered around implying that there is no preferred orientation 
in the configurations and that the system is in isotropic phase. While p > 0.4, the distribution 
shifts to the right indicating onset of a nematic phase. Also, the order parameter increases 
gradually, indicating a continuous phase transition. Since we are not able to determine the 
critical point with very good precision, we did not attempt to study the critical exponents 
characterizing the transition. As there are two competing ordered states near the transition, 
one would expect this transition to be in the Ising universality class. 

The high density transition. - For intermediate densities, the state of the system is well- 
approximated by the fully aligned state in which all the rods are in the same orientation (say 
horizontal). We can determine the entropy of the system approximately, by considering it as 
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Fig. 4 - Normalized distribution of n„ — rih for k = W and L = 120 is shown for different values of 
densities p. 

a perturbation from the fully aligned state. The calculation of the entropy of the fully aligned 
state density p is straight-forward, and reduces to the calculation of a l-dimensional problem. 
The number of ways of putting (pL/k) rods of length A; on a one dimensional line of L sites 
is easily seen to be 



[L(l-p) + lpL]! 
[L(l-p)]!(ipL)! 



(4) 



^From this, it is easily seen that entropy per site in the thermodynamic limit L ^ oo in 
this reference nematic state is given by 



Snemip) = (1 - P + ^) log(l - p + ^) - {1 - p) log(l - p) - ^ log ^ 

For p= 1 — e, this tends to zero as e tends to zero as 

Snemip = 1 — e) = £ log — £ + higher order terms in e 
fee 



(5) 



(6) 



We can expand fl{zh, Zy) in powers of for a fixed value of Zh- We treat Zy as small, and 
only at the end of calculation put z-^ = z^ 



n{zh, z,) = n{zh,m + z^Nh + zt{N'ftl2 + Nh) + 



(7) 



where is the number of sites in the lattice, /i is expressible in terms of the probability 
of finding k consecutive empty sites in the vertical direction in the reference state. Since 
the configurations of different horizontal rods are independent in the reference state, and 
the probability of finding a random site unoccupied is e, we get /i = e*^. Using cumulant 
expansion, it is easy to see that 



5 (p) - S'„em(p) = 2(1 - p)*" + 2:^/2 + higher order terms. 



(8) 
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Fig. 5 - Entropy per site for nematic and disordered states as function of e for fc = 8. 



The value of /i is quite small even for moderately large values of e and small k, Thus for 
e = 1/4 and k = 8, the lowest order term in the perturbation expansion of the fraction of 
vertical fc-mers in the nematic state is only of the order of 10~^. This is much less than about 
1% fraction of vertical bonds seen in Fig. 2. Thus the higher order terms are important. 

The second order term /2 can be expressed in terms of the 2-point correlation functions of 
the unperturbed problem. The detailed form of /2 is somewhat complicated, and is omitted 
here. It can be shown that /2 is of order e'^'^k'', which has the correct order of magnitude. 

However, for e = 0, the fully ordered state is not the most likely state. The number of 
configurations of /c-mers increases exponentially with N of the system as exp[A''S'(p)], where 
^(p = 1) > 0. It is straightforward to see that 



S{p = 1) > 



log 2 

fc2 



(9) 



as there are two ways of covering a k x k square with straight fc-mers, and a lattice of size 
kL X kL can be broken up into such small squares. 

In fact we can get a better bound on 5(1) as follows: break the kL x kL lattice in L strips 
of width k each. Now if Fn is the number of covering an N x k rectangle with /c-mers, we 
easily see that Fn satisfy the recursion relation 



''N-l 



N-k- 



This implies that F^ increases as A where A is the largest root of the equation 

^ ^k-i _^ 1^ 



Form this equation, it is easy to check that for large k, 
A»l 



1 k 

— logf^ -)+ higher order terms in (l/fc). 

k log k 



(10) 



(11) 



(12) 



The entropy per site is then bounded from below by (logA)/fc and this gives the estimate 
5(1) » logfc/fc^ for large k. 

For densities away from the fully packed state an approximate state of the system is 
obtained by first starting with a random configuration of the fully packed lattice, and then 
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remove a fraction of the rods at random. This gives an approximate expression for the entropy 
of the disordered state as 

Sd,s{p = 1 - e) « Sd^s{p = 1) + ^[-eloge - (1 - e) log(l - e)] (13) 

In Fig. [5l we have plotted Sdis{p) and Snem{p) for p near 1. We note that for small e, 
Snemip) increases with e as elog(l/e), but Sdis{p) has a weaker e— dependence, as it varies 
as (e/fc) log(l/e). Clearly the disordered state has higher entropy for p = 1, but for lower 
densities, the nematic state is favored. It is easy to verify that the two curves cross at 
Pc2 ~ 1 — C/k'^, for large k, where C is some constant. Hence we expect a second phase 
transition from the nematic ordered state to the disordered state as p is increased beyond a 
density pc2- 

Unfortunately, we are not able to study this second transition by our Monte Carlo simula- 
tions. While the Monte Carlo algorithm is formally ergodic for all p 7^ 1, the relaxation times 
increases very fast as the density increases and the system gets trapped in some small set of 
nearby states. Our algorithm works well only for not too large values of p ( say p < .85). 

We can only provide a qualitative description of this second phase transition. If we start 
with a fully packed configuration, and remove a single /c-mer, we get a state with k unoccupied 
sites ('monomers'). These monomers occupy consecutive sites in the horizontal or vertical 
direction. If the A:-mers are allowed to diffuse, equivalently the monomers can diffuse, but may 
remain bound together as a molecular bound state of k monomers. However, at higher densities 
of such defects, these monomers become unbound. In the nematic state, the monomers form 
a nearly ideal gas, and the larger entropy of the gas of monomers makes the nematic phase 
preferred thermodynamically. Characterizing this phase transition, in particular to determine 
if it is first order or continuous is an interesting open problem. 
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